Demixing of compact polymers chains in three dimensions 
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We present a Monte Carlo algorithm that provides efficient and unbiased sampling of polymer 
melts consisting of two chains of equal length that jointly visit all the sites of a cubic lattice with 
rod geometry L x L x rL and non-periodic (hard wall) boundary conditions. Using this algorithm 
for chains of length up to 40 000 monomers and aspect ratios 1 < r < 10, we show that in the limit 
of a large lattice the two chains phase separate. This demixing phenomenon is present already for 
r — 1, and becomes more pronounced, albeit not perfect, as r is increased. 

PACS numbers: 82.35.Lr, 05.10.Ln 



INTRODUCTION 

When a single polymer chain with excluded volume 
interactions is brought into a bad solvent, the solvent 
molecules strongly repel those of the chain. To minimize 
the interface between solvent and chain molecules, the 
chain will collapse and expel the solvent molecules from 
its interior. Similarly, a few polymer chains in a very 
large volume of solvent will form a dilute solution of col- 
lapsed chains. In particular, the individual chains will 
not mix [l|. 

Consider now squeezing the chains together so as to 
form a solution in which the monomer density is finite 
(and the solvent density can be taken to zero). We may 
then ask whether different chains will mix or segregate. 
One would expect that the answer to this question might 
depend on how the number of chains scales with the size 
of the system. 

A well-studied case is that of a polymer melt where 
many chains are present. If there are of the order iV^/^ 
chains, occupying a total volume N (i.e., with an average 
monomer density of N^^/^ for each chain) the so-called 
Flory theorem [2[ is supposed to apply. According to 
this, the repulsive hard-core interaction between differ- 
ent chains is completely screened out at large distances. 
This implies that a single chain within the melt exhibits 
ideal, Gaussian behavior: its radius of gyration scales as 
the square root of the chain length. In particular, dif- 
ferent chains should mix (interpenetrate) completely, as 
opposed to a demixing (segregation) scenario where each 
chain would occupy a small lump of the available space 

a- 

It has recently been shown by the Strasbourg group 
in a number of studies [3[ that the Flory theorem is not 
quite correct. In particular, these authors have shown 
that the interplay between chain connectivity and the 
incompressibility of the melt leads to an effective repul- 
sion between chain segments. The bond-bond correlation 
function, supposed to be short-range by the Flory argu- 
ment, turns out to exhibit a long-range algebraic decay. 
A number of deviations from the ideal Flory behavior 



have been pointed out, and shown to be accounted for 
in an improved scaling theory [3| that views the chain 
as a hierarchical arrangement of correlation holes [1| of 
sub-chains on all length scales. 

On the other hand, one might be interested in the sit- 
uation where there is only a finite number of chains in a 
solution of finite monomer density. In this context, the is- 
sue of polymer mixing has recently gained renewed inter- 
est in the area of physical biology. One outstanding ques- 
tion is to understand the organization and segregation 
of chromosomes within simple bacteria, including during 
duplication. In particular, the bacterium Escherichia coli 
with a single circular chromosome in a rod-shaped cell 
has served as a model system. In this situation the ra- 
dius of gyration of a single chain will necessarily scale 
as the cube root of the chain length, and there is no 
good reason why the Flory theory y and the more re- 
cent improvements of it [3| should apply. The question 
whether the chains will mix in that case therefore needs 
to be settled separately, as do the possible consequences 
for chromosome organization. 

The behavior of a dense solution of just two polymer 
chains confined to a bar-shaped geometry, akin to that 
of a rod-shaped cell, was discussed in [j| . A simple pis- 
ton model, combined with a blob-type analysis [l[ (see 
also [5|), was used to point out the mechanisms by which 
demixing might occur. In the limit of large piston pres- 
sure (corresponding to maximal density) the two chains 
were nevertheless found to mix. In contrast with this, 
studies of unentangled ring polymers in concentrated so- 
lutions lD, y concluded that the existence of topological 
barriers would lead to segregation. A molecular dynamics 
study, using a parameter-free bead-spring polymer model 
of a ring polymer melt, corroborated the segregation sce- 
nario and found several structural and dynamical results 
in agreement with experiments on Drosophila and bud- 
ding yeast chromosomes [9[ . The results of [9[ imply that 
a certain number of experimental observations could be 
due to generic polymer effects, and as such could be ac- 
counted for in very simple models. 

In this paper, we revisit the problem of two open chains 



(not rings). In line with the above reductionist point of 
view, we do so by studying numerically a simple and 
very precisely defined parameter-free lattice model of a 
two-chain melt of unit monomer density. This model 
accounts fully for the constraints of maximal monomer 
density, and the self and mutual avoidance of the chains. 

The conformations of a single maximally compact 
(space filling) chain can be modelled on a lattice as a 
Hamiltonian walk (HW). By definition, a HW is a self- 
avoiding walk that visits each lattice site exactly once. 
Upon adding further local interactions, the HW model 
has been proposed as a description of protein melting [g 
(with bending rigidity), or of protein folding [10[ (with 
suitable interactions among amino acids). The simplest 
example of the latter is the so-called HP model in which 
only two types of amino acids (Hydrophobic or Polar) are 
taken into account. 

To similarly model a melt of a finite number of chains, 
we can use the model of Hamiltonian chains (HC) [ll| in 
which M self- and mutually avoiding walks jointly visit 
all lattice sites. By definition, each chain has length at 
least one (its two end points cannot coincide). In what 
follows we shall always take the underlying lattice to be 
simple cubic (or square, when we occasionally compare 
with the two-dimensional case). It seems natural to as- 
sume that the HC model correctly describes an M-chain 
polymer melt in the large-scale limit, where lattice details 
should be irrelevant. 

The question is then whether the HC model validates, 
for the case of a melt consisting of only a few chains, the 
mixing scenario envisaged for the two-chain case in [j]. 
In this paper we study in detail the simplest possible case 
of M = 2 identical chains of equal length, by large-scale 
numerical simulations. We show that, contrary to the 
above expectations, the chains actually phase separate. 
For an i X L X (rL) lattice, with large L and fixed aspect 
ratio r, the effect turns out to be quite subtle (but mea- 
surable) for r = 1, becoming more pronounced for larger 
values of r. 



ALGORITHM 

Properties of self-avoiding walks are routinely stud- 
ied by importance sampling schemes of the Monte Carlo 
(MC) type, one well-known example being the pivot algo- 
rithm J12l |. Unfortunately, almost all known algorithms 
become useless in the fully-packed limit, as the accep- 
tance ratio for any proposed move tends to zero. For the 
HW case (i.e., the HC model with M = 1 chain) we have 
recently shown that an MC move proposed by Mansfield 
[13| can be turned into an algorithm that is both efficient 
(with dynamical exponent z ~ 0) and provides unbiased 
sampling (each conformation is generated with uniform 
probability) [l3|. 

The working principle of this M — 1 algorithm is 
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FIG. 1: (Color online.) Move between two of the 2 480 304 
possible Hamiltonian walks on a 3 x 3 x 3 cubic lattice. 



shown in Fig. [T] Choose an end point x of the HW at 
random, and select one of its adjacent empty links (xy). 
Then, among the occupied links adjacent to y, there is 
exactly one which would form part of a loop if (xy) were 
added to the HW; call it (yz) . The basic move is then to 
change the HW by adding (xy) and removing (yz). 

We now state the modifications necessary to deal with 
the case of M = 2. For the moment we assume that the 
point x is in the bulk of the system; the case when it is 
on the boundary will be discussed later. 

Let the configuration C consist of two chains, and let 
x be one of the four end points chosen at random. Once 
again, choose with uniform probability one of its adjacent 
empty links (xy). As before, the idea will be to modify C 
by adding (xy) and removing one other link. But before 
adding the link (xy), the point y may be adjacent on 
one or two occupied links. The former happens if y is 
another end point, either of the same chain as x or of the 
other chain; in that case the choice of the link (yz) to be 
removed is unambiguous. In the latter case (i.e., when y 
is adjacent on two occupied links) one must distinguish 
two different situations. If y is on the same chain as x, 
one proceeds as in the M = 1 case, making the unique 
choice of z that avoids the formation of a loop. And 
if x and y are on different chains, one choses randomly 
between the two possibilities for z. 

It remains to discuss two subtleties. First, when x 
and y are end points of two different chains, the move 
will consist in letting the first (resp. second) chain grow 
(resp. retract) by one monomer. In particular, it might 
happen that z is the other end point of the second chain. 
The move would then consist in letting the second chain 
shrink to zero, but since this is not an allowed HC config- 
uration, we remedy the situation by leaving C unchanged. 
Second, when x is at the boundary of the system, it might 
happen that y is chosen outside the lattice. In that case 
as well the move is rejected, i.e., C is left unchanged. 



UNBIASED SAMPLING 

With all these rules, and considering that an MC move 
takes unit time regardless of whether it leads to a change 



in C or not, it follows from a careful analysis that the 
rule of detailed balance is satisfied. As an independent 
numerical check, we ran the simulation on a 2 x 2 x 3 cubic 
lattice until each of the 4 204 possible configurations [11| 
had been generated ~ 5 • 10^ times. We then verified 
that the fluctuations in the number of occurencies of each 
configuration was compatible with counting statistics. 

Proving that the MC move is also ergodic, i.e., that 
any desired configuration may be reached from an initial 
one in a finite number of steps, unfortunately turns out 
to be much harder (a proof does not even exist for the 
M = f chain case 14|). 

As in 1J| , we therefore resort to explicit checks of er- 
godicity for small 2D and 3D systems. The case of a 
2x2 square lattice is special, since no MC move con- 
nects the two possible configurations. But for all other 
Lx X Ly lattices with 2 < Lx,Ly < 6, we have checked 
that the number of different configurations generated by 
the MC move, starting from an initial reference config- 
uration, coincides precisely with the exact enumeration 
results of [ll|. In particular, precisely 17 570 172 config- 
urations are found on the 6x6 lattice. We have similarly 
checked ergodicity on 2 x 2 x 2, and 2x2x3, and 2x3x3 
cubic lattices. 

Assuming that ergodicity holds in general (except for 
the trivial 2x2 square lattice) , we conclude that the MC 
process converges towards the equilibrium distribution, 
i.e., that all HC configurations arc sampled with uniform 
probability. 



PROPERTIES 

In the case where the point y is an internal point on 
the same chain as the end point x, the identification of 
the correct edge (yz) to be removed necessitates tracing 
out the loop formed when adding xy. So at worst, one 
MC move may require a time '^ N. Since each move 
changes at most two links, it is appropriate to evalu- 
ate the autocorrelation time r (defined in terms of the 
link overlap with the initial configuration) in units of the 
number of MC moves per site. Just as in the M = 1 case 
[l4J we find that in those units r ~ L^ is of order unity, 
independently of the system size. We conclude that the 
dynamical exponent 2 « 0, proving our claim that the 
algorithm is efficient. 

Note that our algorithm provides unbiased sampling of 
configurations of two chains of lengths A^i and 7V2 , where 
only A^ = A^i -I- A^2 is fixed. To pursue our goal — the 
study whether two chains of equal length mix or not — 
we must impose the constraint A^i = A''2 in an efficient 
way. Trial runs for L = 16 and aspect ratios r = 1,2 
show that the probability distribution p{p) of the length 
ratio p = Ni/N is very close to being uniform: we have 
in fact \p{p) - 1| < 0.01 for ah p S [0.02, 0.98]. It follows 
that a configuration with p = 5 '^'^^^ occur on average 




FIG. 2: (Color online.) A randomly chosen configuration of 
two chains of equal length on an 8 x 8 x 16 cubic lattice. 



once for each unit of t. 

It is thus tempting to collect a data point whenever 
/9 = i. There is however an important caveat. Suppose 
that the end points of chain 1 are deeply immersed in 
chain 1, and same for chain 2, and that the two chains 
happen to be of the same length. Then all MC moves 
will be such that an end point x of some chain attacks 
a point y on the same chain, and never on the other, 
thus conserving the length of either chain. If data were 
collected whenever p = ^, a, lot of measurements would 
be made which were spaced by just one MC move, and 
not one unit of r. This would strongly bias the results. 
A better measuring protocol is to collect a data point 
when p = 2 '^"■^ require two data points to be spaced 
by a least one autocorrelation time. With this protocol, 
each data point gives an independent configuration with 
p = ^ , taken uniformly within the ensemble of chains of 
equal length. 

To avoid transient effects due to the choice of the initial 
state, we discard in each run the first 100 MC moves per 
site, before starting the data collection. Fig. [5] shows a 
measurement on a rather small system obtained using 
this protocol. To the eye it certainly looks like the two 
chains interpenetrate, but obtaining a valid conclusion 
will require analysing the system more carefully. 



SCALING THEORY 

According to the standard scaling theory of polymers 
m, the number of HW of length A^ is - fi'^ N^-^ for 
N ^ 1. Here /i is a lattice-dependent connective con- 
stant, and 7 a universal conformational critical exponent. 
Alternatively, 7 is determined by P{x) ~ x^"*", where 
P{x) is the probability distribution of the end-to-end dis- 
tance in the asymptotic regime a; <C 1. Using this, the 
conformational exponent for HW has been determined 
numerically as 7 = 0.94 ± 0.02 



14| . More recently, the 



consistent value 7 = 0.95 ± 0.02 has been reported [li 

The value of 7 can be used in various ways to ar- 
gue whether M = 2 chains will phase separate or not. 
One simple argument goes as follows. Suppose the two 



chains live in a volume V. In the demixing scenario, 
each chain would occupy a volume V/2, and the number 
of configurations would scale like that of two indepen- 
dent HW, viz. [^^/2(y/2)T-i]^ - //^y2(7-i), if^ on the 
other hand, the two chains mixed, one could make the 
assumption (A) that one end point of either chain would 
be 'close' and hence could be connected so as to form a 
single HW in volume V. The number of configurations 
would then be ii^V''~^. Comparing these, we conclude 
that entropy should drive the chains to mix if and only 
if 7 < 1. The above numerical determinations of 7 are 
then in favour of the mixing scenario. 

The applicability of the above argument clearly hinges 
on whether the assumption (A) can be considered cor- 
rect. Some partial justification can be found in |ll| in 
which the interactions of a single HW with the surface 
of the system was carefully analyzed. In particular, it 
was found that the end points possess a clear tendency 
to avoid lattice sites of high coordination number, i.e., 
they are attracted towards the surface. Since there are 
fewer surface sites than bulk sites, this implies an en- 
hanced probability for end points to be close than would 
have been the case if they were uniformly distributed 
throughout the volume. Moreover, the fact that 7 < 1 
means that end points attract one another entropically 
also in the bulk. 

But despite of such arguments, the assumption (A) 
is certainly a rather weak spot in the above argument: 
the existence of an attraction between end points does 
of course not mean that they are separated by a single 
lattice spacing. So the above reasoning cannot be con- 
sidered as conclusive. The same is true for other variant 
scaling arguments that we have tried out. We therefore 
turn to our numerical results to determine whether the 
two chains will mix. 



NUMERICS 

We have used our MC algorithm to generate an exten- 
sive number of independent configurations for the two- 
chain problem on L x L x {rL) cubic lattices with L = 16 
or 1/ = 20 and various aspect ratios 1 < r < 10. In 
the largest simulation (L = 20 and r — 10) each of the 
chains had length ^L^ = 40 000 monomers and up to 
1.1-10^ statistically independent configurations were gen- 
erated. The boundary conditions are non-periodic (hard 
wall confinement) in all three lattice directions, as shown 
in Fig. m 

In a first series of runs, we produced A/i ~ 3 • 10^ 
configurations for each r. For each of these, the standard 
- (z)2 of the ^-coordinates of the 
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monomers — rescaled so that z S [0, 1] — was computed 
separately for each chain, and an average value (cz) was 
extracted from the 27Vi measures. The results for L = 16 
systems arc given in the left part of Table ID They can be 
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0.28231 ± 0.00005 
0.27675 ± 0.00008 
0.27094 ± 0.00011 
0.26527 ± 0.00019 
0.26002 ± 0.00015 



0.3632 ± 0.0005 
0.3580 ± 0.0005 
0.3560 ± 0.0005 
0.3550 ± 0.0003 
0.3546 ± 0.0003 



TABLE I: For various aspect ratios r, we show: (1) the aver- 
aged standard deviations Ip-z) of the z-coordinate of a chain 
in tlie fully-packed two-chain problem, and (2) the monomer 
fiicacitv 7/iu. that mil make a. sine'le chain fill out one half of 




FIG. 3: (Color online.) Probability distribution p{/S.z) of the 
difference of center of mass coordinates Az = \z\ — Z2\ of the 
two chains, for L = 16 lattices with various aspect ratios r. 
The dashed curve is for L — 2Q and r = 10. 



compared to the result a^ = -K= ~ 0.28868 which would 
be obtained in the case of full mixing, where the mass of 
each chain is uniformly distributed in the z-direction. 

The results for {uz) indicate that the two chains tend 
to demix, and that this tendency grows with increasing 
r. Of course, one is far away from the value ^crl which 
would result in the hypothetical situation of complete 
demixing, where one chain goes completely to the left of 
the box, and the other to the right. Note also that the 
difference from a^ is clearly measurable — albeit tiny — for 
r = \ where the three coordinates are equivalent. 

To give more evidence for the demixing, and charac- 
terise it more precisely, we turn to a second series of runs. 
In this series we construct the probability density p{l^z) 
of the quantity Az = |zi — Z2I, where Zi = (z) is the 
z-coordinate of the center of mass of the monomers in 
chain i (for i ~ 1,2). Obviously, constructing a whole 
probability distribution calls for better statistics, so in 
this series we produced A/'2 — 1.1 • 10^ configurations for 
each r. The results for p{Az) are shown in Fig. [S] To 
eliminate fluctuations arising from the lattice discretiza- 
tion, the values of Az have been arranged into 50 'run- 
ning' bins; the statistical error bars can be judged from 
the remaining small ripples on the curves. 
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FIG. 4: (Color online.) Mass density pii{z) of one of the 
chains, measured relative to its center of mass coordinate zi. 
The data are for L — 20 lattices with various aspect ratios r. 
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FIG. 5: (Color online.) Mass density pi2{z) of chain 1, mea- 
sured relative to the center of mass coordinate 22 of the other 
chain. The data are for L — 20 lattices with various aspect 
ratios r. 



In a third scries of runs we have focused on the mass 
density Pij{z) of chain i, measured relative to the center 
of mass Zj of chain j (for i = 1,2). The normalization 
is chosen so that J pij{z)dz = 1. Demixing can be de- 
tected via a different behavior for the density of a chain 
in its own center of mass (pn — P22 by symmetry) and 
for the density of a chain in the center of mass of the 
other chain (pi2 = P2i)- For these measurements we have 
taken L = 20 and produced up to A/3 ~ 3 • 10'' configu- 
rations for each r. Instead of using bins, we have simply 
rounded the center of mass coordinates as expressed in 
lattice spacings (i.e., rLzj) to the nearest integer. The 
results for pii{z) are given in Fig. HI and those for ^12(2:) 
in Fig. [5] Both functions should of course be symmetric 
upon reversing the coordinate axis, Pij{z) = pij{—z), and 



the error bars can be judged from the slight deviations 
from this symmetry. It is seen from the figures that when 
r increases, pii{z) tends to become a more narrow dis- 
tribution around the origin, whereas pi2{z) spreads out 
and develops a depletion near the origin. This is again 
clear evidence of demixing. 



DISCUSSION 

The data in Fig.[3]give clear evidence for demixing as r 
increases. For small r we see a rather broad distribution 
p(Az), localized around the origin. At r ~ 7, the distri- 
bution becomes almost uniform, and for higher r a peak 
near Az — i, the maximal possible value, develops on 
a close-to-uniform background. To make sure that this 
conclusion is not a finite-size artefact, we recomputed 
the r = 10 curve for a larger size L = 20. Note that this 
computation required several years of CPU time. The re- 
sult, shown as a dashed curve in Fig. [U does not deviate 
substantially from the L = 16 result. 

It is useful to compare this situation to that of a sin- 
gle non-fully packed chain that is constrained to take up 
precisely one half of the available volume. To study this, 
we employ the non-fully packed version of the M = 1 al- 
gorithm, discussed very briefly at the end of [1J|. To set 
the peak of the monomer concentration p = N/V pre- 
cisely at p = i , one needs to tune the Boltzmann weight 
w of a monomer. Still for L = 16, we find that this is 
obtained for w ^ w^,, where values of Wi, for different 
r are given in the right part of Table ID This allows us 
to take a series of data for which each configuration has 
p = ^ precisely. Based on a set of A/s ~ lO'^ data points, 
we then construct the probability distribution p{zi) of 
the center of mass zi of the chain. For convenience, we 
here subtract i from the z-coordinates so as to set the 
origin in the middle of the box; the support oi p{zi) is 
then zi S [-"I:, -j]- The results for p{zi), using this time 
200 bins, are given in Fig. [SI 

It is of course completely predictable that these distri- 
butions arc Gaussian; we find indeed p{zi) ~ aexpf — 
6(zi)^) with b = (0.91 ±0.01) x r. But more importantly, 
comparing Figs. [3]-[6l one can appreciate that even for 
r = 1 the distribution p(Az) is indeed very broad. This 
corroborates our previous claims that even for r = 1 the 
two-chain problem indeed exhibits demixing. 

It is useful to compare the data in Figs. HHS] to the 
hypothetical situations of complete mixing (with both 
chains uniformly distributed in the z-direction) and of 
complete demixing (with one chain occupying the left half 
of the box, < z < i , and the other chain the right half, 
2 < z < 1). In these two extremal situations, the mass 
distributions pii{z) and pi2{z) would be those shown in 
Fig. [Zl The data of Figs. HHS] bear some resemblance to 
a crossover (with increasing r) between these extremes. 
A closer examination of the data however indicates that 
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FIG. 6: (Color online.) Probability distribution p{zi) of the 
centre of mass coordinate 21 of a single chain, constrained to 
occupy precisely one half of the box, for different aspect ratios 
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CONCLUSION 

We have studied a problem of two compact polymer 
chains of equal length in a three-dimensional box of var- 
ious aspect ratios r. We have devised an efficient Monte 
Carlo method for sampling all configurations of this prob- 
lem with uniform probability. Naive scaling theory, com- 
bined with existing evaluations of the conformational 
exponent 7 for a Hamiltonian walk, indicates that the 
two chains should mix completely. This same conclusion 
would be suggested by the Flory theorem [6| for a poly- 
mer melt of many chains, or from a blob-type analysis for 
a simple piston model of just two chains [j|. However, our 
numerical simulations clearly show that the two chains 
do in fact phase separate, even for r = 1, and that this 
demixing becomes more pronounced, albeit not perfect, 
in the limit of large r. 

Clearly, further work is needed to see whether this con- 
clusion can be reconciled with soirre iirrproved scaling the- 
ory, specifically adapted to a melt consisting of just a few 
chains. The demixing scenario obtained in this paper is 
sonrewhat subtle, since it is not due to the additive terms 
of the free energy. Further study is therefore required 
to assess whether the effect is stable towards introduc- 
ing more parameters into the model (non-identical chain 
lengths, specific chemical interactions at the chain ends, 
off-lattice modelling, etc.). It would also be interesting 
to see whether demixing can be realized experimentally. 
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FIG. 7: (Color online.) The appearance of the distributions 
pii{z) (left panels, cf. Fig. [4| and pi2{z) (right panels, cf. 
Fig. [5]) in the two extremal and hypothetical situations of 
complete mixing (top panels, red or light gray curves) and 
complete demixing (bottom panels, blue or dark gray curves). 



for r 3> 1 the distributions pij (z) will not coincide with 
those of the complete demixing scenario. (For instance, 
the tails in Fig. |4] extending out to 2; = ±| seem to be 
robust upon increasing r.) Note also that in terms of the 
distribution p{Az) of Fig. [3] complete demixing would 
correspond to a Dirac delta function S{z — ^), whereas 
the actual data are indicative of a close-to-uniform dis- 
tribution on [0, 2], plus possible a finite peak near z = ^. 
In other words, the numerics strongly supports the idea 
that the physics in the limit limr_>.oo i^i^L->oo (with the 
thermodynamic limit being taken first) is non-trivial and 
corresponds to some kind of "partial demixing" . 
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